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We present two complementary models to study stationary nonlinear solutions in one-dimensional 
plasmonic slot waveguides made of a finite-thickness nonlinear dielectric core surrounded by metal 
regions. The considered nonlinearity is of focusing Kerr type. In the first model, it is assumed 
that the nonlinear term depends only on the transverse component of the electric held and that the 
nonlinear refractive index change is small compared to the linear part of the refractive index. This 
hrst model allows us to describe analytically the held prohles in the whole waveguide using Jacobi 
elliptic special functions. It also provides a closed analytical formula for the nonlinear dispersion 
relation. In the second model, the full dependency of the Kerr nonlinearity on the electric held 
components is taken into account and no assumption is required on the amplitude of the nonlinear 
term. The disadvantage of this approach is that the held prohles must be computed numerically. 
Nevertheless analytical constraints are obtained to reduce the parameter space where the solutions 
of the nonlinear dispersion relations are sought. 

PACS numbers: 42.65.Wi, 42.65.Tg, 42.65.Hw, 73.20.Mf 


I. INTRODUCTION 

Studies of stationary nonlinear waves possessing the 
properties of both plasmons and solitons started in the 
early 80s when this type of waves was described by Agra¬ 
novich et al. [T]. Two main types of structures supporting 
plasmon-soliton waves were studied. The first type con¬ 
tains one or two semi-infinite nonlinear media and was 
extensively studied in Refs. HHE]. Transverse electric 
(TE) and transverse magnetic (TM) polarized waves were 
investigated in configurations built of semi-infinite non¬ 
linear layers, a metal layer, and, in some cases, additional 
dielectric layers. The second type of structures contains a 
finite-size nonlinear dielectric layer sandwiched between 
two semi-infinite metal layers. This type of structures 
will be called here nonlinear slot waveguide (NSW). Stud¬ 
ies of structure with a nonlinear dielectric core started in 
the early 80s from fully dielectric structures [HHIS]. The 
solutions of Maxwell’s equations in a finite size nonlin¬ 
ear dielectric were given in terms of Jacobi elliptic func¬ 
tions |24] . A symmetry breaking bifurcation was pre¬ 
dicted for the fundamental symmetric mode giving birth 
to an asymmetric mode [m US]. Various methods to 
study fully dielectric waveguides with both nonlinear core 
and cladding were proposed 

In 2007, Feigenbaum and Orenstein [35] made the first 
attempt to study waveguides with a nonlinear core sur¬ 
rounded by a metal cladding (NSWs) instead of a dielec¬ 
tric one. Their method describes sub-wavelength confine¬ 
ment of light in two-dimensional plasmon-soliton beams 
propagating in NSWs. Such a strong confinement is en¬ 
sured by a linear plasmon profile in the transverse direc- 
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tion and by the self-focusing effect in the lateral direc¬ 
tion. Recently, numerical [M] and semi-analytical [531133] 
methods have been developed to study up to the three 
first nonlinear modes in NSWs. Higher order modes were 
also reported in NSWs [57] . 

The NSW is an interesting and promising configuration 
for two reasons: (i) from the practical point of view, it is 
easier to fabricate high quality thin nonlinear films than 
bulky layers, like those needed in the configurations with 
semi-infinite nonlinear medium and (ii) it has numerous 
potential applications. Devices based on the NSW con¬ 
figuration can be used for phase matching in higher har¬ 
monic generation processes [38] and for nonlinear plas¬ 
monic couplers [39j|40]. Nonlinear switching [41] was the¬ 
oretically predicted in NSW-based structures that is sim¬ 
ilar to the nonlinear switching in graphene couplers [42] . 
Tapered NSWs might be used for nanofocusing and loss 
compensation in order to enhance nonlinear effects [43] . 

The field of NSWs is relatively young and there is not 
a lot of works describing the properties of these struc¬ 
tures. In this article, we build two new and complemen¬ 
tary models that allow us to study efficiently and accu¬ 
rately the NSW configurations. These models have been 
briefly introduced in Ref. m and here we present their 
detailed derivations. 


In Sec. jlTJ a general description of the problem studied 
here is given. Section III presents the theoretical deriva¬ 
tion of our two models. In Sec. |IV[ the validity of the 
model formulations is verified in the limiting case of a 
semi-infinite medium. These models are used to thor¬ 
oughly study the propagation of plasmon-soliton waves 
in NSWs in the following article [44] where analytical and 
numerical stability analysis results are also provided. 
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II. PROBLEM STATEMENT 

This article presents a derivation of the dispersion re¬ 
lations for the stationary TM polarized waves propagat¬ 
ing in one-dimensional NSWs depicted schematically in 
Fig. [3 We propose two models to study the light prop¬ 
agation in such structures. The hrst model is based 
on the approach proposed for fully dielectric structures 
in [ninD] but is extended to structures containing met¬ 
als. Due to the presence of the metal cladding, the field 
continuity conditions result in a new type of held prohles 
and very rich dispersion relations. This approach uses an 
approximated treatment of the nonlinearity in the Kerr 
medium which allows us to write and solve a single non¬ 
linear wave equation in the hnite-size nonlinear medium. 
Using the held continuity conditions at the core interfaces 
located at x = 0 and x = d the analytical formulas for the 
dispersion relations and for the held prohles are obtained 
in terms of Jacobi elliptic functions [23]. Therefore, this 
model will be called Jacobi elliptic function based model 
(JEM). 

The second model is based on the approaches from 
Refs. [21 Hn] for a single interface between a nonlinear 
dielectric and a metal, and therefore it is named the 
interface model (IM). This model uses a more realistic 
treatment of the nonlinear Kerr effect than the JEM {i.e. 
Kerr nonlinearity depends on all the components of the 
electric held). It allows us to obtain separate dispersion 
relations on the two interfaces of the NSW in analytical 
forms. Comparing the dispersion equations for the left 
interface and for the right interface, results in an analyt¬ 
ical condition that reduces the parameter space in which 
the solutions of Maxwell’s equations in NSWs are sought. 
The solutions are found by the numerical integration of 
Maxwell’s equation in the core which allows then to re¬ 
late the two interfaces. Maxwell’s equations in the core 
are solved using the shooting method [45]. If the result 
of integration is consistent with the previously assumed 
values of the held and its derivatives at the slot inter¬ 
faces then the corresponding j3 is accepted as a genuine 
solution in our problem. 

The stationary solutions in our one-dimensional geom¬ 
etry are sought in the form of monochromatic harmonic 
waves: 
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where E = and H = [0,Hy,0]. We write 

the imaginary unit i in front of the z-component of the 
electric held so that all the quantities E^,Ez, and Hy 
are real. The propagation direction is chosen to be z and 
w denotes the angular frequency of the wave, fcg = ujjc 
denotes the wavenumber in vacuum, c denotes the speed 
of light in vacuum, and jd denotes the effective index of 
the propagating wave. The structure is invariant along 
the y direction and therefore, it is assumed that the held 
prohles are invariant along the y coordinate. 


III. MODEL DERIVATION 
A. Maxwell’s equations 


The derivation of our models starts from the gen¬ 
eral form of the Maxwell’s equations in case of nonmag¬ 
netic materials (relative permeability y, = 1) without free 
charges {pf = 0) and free currents (J/ = 0) [35]. Eor the 
harmonic monochromatic TM waves described by Eq. Q 
and for an isotropic relative permittivity e the Maxwell’s 
equations read 


kol3Ex - = ujpoHy, 


E^ = 


-H„ 


egec 


E, = 


1 dHy 


€oeoj dx 

The nonlinearity studied here is of the Kerr type 


(2a) 

(2b) 

(2c) 


e = ei{x) + eni(x), 


(3) 


where ei denotes the linear, real part of the permittivity, 
Eni = 0^(2’) IE ( 2 :) P denotes the nonlinear part of the per¬ 
mittivity limited to the isotropic optical Kerr effect that 
depends on the electric held intensity, and a{x) denotes 
the function that takes values of the nonlinear parame¬ 
ters associated with different layers (in linear materials 
it is null). The imaginary part of the permittivity is ne¬ 
glected in the modal studies presented here. 

The formulation of the Kerr effect used in the follow¬ 
ing of this study can be found in the majority of the 
works on the nonlinear waveguides and nonlinear sur¬ 
face waves [IHzl IIM Iliffl IM® lIZHZl] . It describes 
sufficiently well the nonlinear effects studied here. To de¬ 
scribe the physics of our system we do not need to use 
the more complex form of the Kerr nonlinear term which 
was described in Refs. [751176]. 


B. Jacobi elliptic function based model 

1. Nonlinear field profiles 


FIG. 1. Geometry of the plasmonic NSW with the parameters We start the presentation of the models for NSWs with 
of the structure. the approach that uses strong assumptions on the form 
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of the nonlinear Kerr term but it provides the dispersion 
relations and the field profiles in analytical forms. This 
model provides more insight and understanding of the 
nature of the problem of finding stationary solutions in 
NSWs than the second, more numerical model. First, we 
will solve the nonlinear wave equation inside the waveg¬ 
uide core and find the nonlinear field profiles. Knowing 
the field profiles, we will be able to derive the dispersion 
relations for the NSW using the continuity conditions on 
the two nonlinear core interfaces. 

In the frame of the Jacobi elliptic function based model 
(JEM), the Kerr-type nonlinearity is not treated in an ex¬ 
act manner. We assume that the nonlinear response of 
the material is isotropic and depends only on the trans¬ 
verse component of the electric field in the following 
way miiaisT]: 

e{x) = ei{x) + a{x)E'l{x). (4) 

Functions ei{x) and a{x) are step-wise functions which 
take the values indicated in Table |l] depending on the 
layer (see Fig.for layer number). 


Layer 

Abscissa 

ei{x) 

a{x) 

1 

X < 0 

£1 = (- 1,1 

0 

2 

0 < X < d 

^ 1,2 

(2) 

eocei^ 2 n 2 = 0^2 

3 

X > d 

£3 = £i,3 

0 


TABLE 1. Values of the functions ei{x) and q;(x' 

the properties of the materials in different layers 

( 2 ) 

linear parameter n 2 in layer 2 is denoted by ■ 

describing 
The non- 

The derivation of the JEM starts from Maxwell’s equa¬ 
tions [Eqs. (^]. These equations are combined together 
with Eq. (E) and with the help of the approximations 
about small nonlinear permittivity change we obtain the 
nonlinear wave equation: 

A2 TT 

-^-klq^ix)Hy + kla{x)H^y=0 

(5) 

where 


q2(x) = /32 _ ei{x), 

(6) 

a(x) = /3^a(x)/leo€i(x)cf. 

(7) 


A detailed derivation of Eq. ([^ is presented in Refs. [Ml 
[77]. We use the first integral treatment approach ISllM] 
and integrate Eq. ([^ with respect to x. The result reads 

=Co- (8) 

The left-hand side of this equation gives us a formula for a 
quantity that is conserved along the transverse profile of 
the core of our one-dimensional nonlinear waveguide. Re¬ 
gardless of at which position x we calculate it, the result 
will always be equal to the integration constant cq . Eor 


the structures with a semi-infinite nonlinear medium m, 
the integration constant was set to zero, because both the 
magnetic field Ely and its derivative dEly/dx tend to zero 
as X —>■ ±oo. Therefore, in semi-infinite cladding layers, 
we must set the integration constant cq = 0. Addition¬ 
ally, in these linear layers a(x) is equal to zero. Thus, in 
the cladding, Eq. ([^ reduces to a standard linear wave 
equation whose solutions are given by: 

Hi = for — oo < x < 0, (9a) 

H 3 = for d < X < -foo, (9b) 


where only the appropriate exponential solutions are con¬ 
sidered. Here the magnetic field amplitudes at the in¬ 
terfaces X = 0 and x = d are denoted by Hq and Hd, 
respectively and qk denotes a constant value of the q{x) 
function in fc-th layer (for k G {1,2,3}). As Hy is the 
only component of the magnetic field, in the following 
derivation we omit the y subscript and instead we use 
a subscript that enumerates the layer in which the field 
profile is defined (see Eig.[^. 

The integration constant cq can be expressed as a func¬ 
tion of the magnetic field amplitude Hq using the conti¬ 
nuity conditions for the tangential electromagnetic field 
components {Hy, E^) at x = 0 in Eq. 


Co — kg 


{ei,2/eifql - qi + a2Hl/2 


H 


0 5 


( 10 ) 


where 02 denotes a constant value of a(x) function in the 
nonlinear core (layer 2) and where, based on the assump¬ 
tion that the nonlinear permittivity change is small, we 
have substituted e2|a;=o+ by ei ^2 in the fraction numera¬ 
tor. Similar expression can be obtained for the second 
interface x = d. Equations ® and (101 allow us to find 


the sign of the integration constant cq for certain types of 
solutions. The solutions can be later classified according 
to the sign of cq . Looking at Eq. ® we notice that for Hy 
field profiles that cross zero, at the point where Hy = 0, 
the only nonzero term on the left-hand side of this equa¬ 
tion is (diJy/dx)^, which is strictly positive. Therefore, 
for this type of solutions, cq can only be positive. Erom 
Eq. ( [l0| ) we notice that for negative q\ the value of the 
integration constant cq is also positive. 

We are searching for guided waves in three-layer struc¬ 
tures. Looking at Eqs. , we notice that the condition 
for the waves to be localized in the waveguide core is 
satisfied when both qi and q^ are real and positive quan¬ 
tities. In order to satisfy this condition, we will look only 
for the solutions with /3 > max{ei,e 3 } [see the definition 
of q{x) and q^ given by Eq. ([^]. The quantity q 2 can be 
either real or imaginary leading to positive or negative 
values of q^. 

In order to find the solutions of the nonlinear wave 
equation [Eq. ([^] in the nonlinear core, we rewrite its 
first integral [Eq. ®] in the form: 


dH2 


^co + klqlH^ - kl^H^ 


= ± dx. 


( 11 ) 
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which is then transformed to 

— / ^ = = ± \ — ax, 

+ AQHl - V A 

using the reduced parameters Q and A: 

Q = klql, 

A = {kla2l2)-\ 


( 12 ) 


(13a) 

(13b) 


In this study, we deal only with the focusing Kerr-type 
dielectrics, therefore A is always positive. Parameter Q 
can be either positive or negative depending on the sign 
of ql- 


Solutions of Eq. (12) take different forms depending on 


the sign of parameters cq. We will solve this equation in 
two cases depending on the sign of Cq. The details of this 
derivation are provided in Appendix [X} 

I. The case cq > 0 

At first, we consider the case where cq > 0. The solu¬ 
tion of the nonlinear wave equation [Eq. ([^] for this case 
is given by 


H 2 {x) = (5cn I(x — a;o)|m| 


(14) 


where cn(u|m) is a Jacobi elliptic function with an argu¬ 
ment u and a parameter m [23] , and 


m = J^/s, 
s = 7 ^ -I- 

7^ = (VA2Q2+471co- AQ)/2, 

= (VA 2 Q 2 + 471 co + Ag)/2, and 
Xo = A/sctT^ [i72(0)/J|m] . 


(15a) 

(15b) 

(15c) 

(15d) 

(15e) 


Here cn ^(u\m) is the inverse of the Jacobi elliptic func¬ 
tion Cn(M|TO). 

II The case Co < 0 

Eor Co < 0, the solution of the nonlinear wave equation 
[Eq. ([^] is given by 

H2ix) = 7 dn|\/ 72 /A (x - xo)|to| , (16) 

where dn(u|m) is a Jacobi elliptic function [24] and 



(17a) 

(AQ + v/A2Q2_47i|co|)/2, 

(17b) 

(AQ-^/A2Q2_47i|co|)/2, 

(17c) 

— \/A/ 72 dn“-^ [id2(0)/7 m] . 

(17d) 


Here dn~^(u|m) is the inverse of the Jacobi elliptic func¬ 
tion dn(rt|m). Knowing the analytical form of the solu¬ 
tions in each layer of the structure [Eqs. and (14) or 
(16)] we can derive the nonlinear dispersion relations for 
the cases of positive and negative cq values. 

2. Nonlinear dispersion relations 


I. The case co > 0 

Firstly, we derive the dispersion relation for NSWs 
when Co > 0. Using the analytical formula for the field 
profile in the nonlinear core [Eq. JT^j, the field profiles in 
metal claddings given by Eqs. and Maxwell’s equa¬ 
tions [Eqs. we write the continuity conditions for 
the tangential electromagnetic field components Hy and 
Ez at the interfaces between the nonlinear core and the 
metal cladding x = d {the continuity conditions at x = 0 
was already used to evaluate the integration constant [see 
Eq. ^]} 

1. The continuity condition for the magnetic field 
component 


iJ. 


2 \x—d 


- =H. 


3\x—d+ 


yields 


Jen \/s/A [d — xq) 


= Ea¬ 


rn 


(19) 


2. The continuity condition for the tangential electric 
field component 


E 


z,2\x—d~ 


= E. 


'z,3\x—d'^ 


( 20 ) 


transformed with the use of Eq. (2c), the formu¬ 


las for the Jacobi elliptic function derivatives, and 
the symmetry properties of Jacobi elliptic functions 
gives 


Je 3 V^/(fcog 3 Q, 2 VA) sn \fsjA (d-xo) I m 


dn 


\/s/A (d- Xo) I m 


= Hd, 


( 21 ) 


where, based on the assumption that the nonlinear 
permittivity change is small, we substituted e 2 \x=o+ 
by £(,2 in the denominator on the left-hand side. 


Comparing the two expressions for EI^ given by Eqs. (19) 


and (21), we obtain the nonlinear dispersion relation in 


its final form for the case of cq > 0 


fco 93 £i. 2 '\/Acn sjA (d — xq) I m = e 3 fcog 3 ei, 2 -\/isn a/ sjA (d — xq) I vn dn a/ sjA (d — xq) 


( 22 ) 
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II. The case Co < 0 

Here we derive the dispersion relation for NSWs when 
Co < 0. The method used here is exactly the same as in 
the previous case (for cq > 0). Using the analytical for¬ 
mula for the field profile in the nonlinear core [Eq. (16)] 
the field profiles in metal claddings given by Eqs. ^1, 
and Maxwell’s equations [Eqs. (1^], we write the conti¬ 


nuity conditions for the tangential electromagnetic field 
components Hy and at the interfaces between the 
nonlinear core and the metal cladding x = d The two 
resulting expressions for Hd are compared and give the 
nonlinear dispersion relation in its final form for the case 

Co < 0 


koq3ei,2'/Adn ^/jyA{d - Xq) 


/r 


O / A / 1 






Equations (221 and (23) build the full dispersion relation 
for the NSW. In order to obtain the dispersion diagram 
for a fixed structure and wavelength (ci, ei^ 2 , ^ d, 

A) we scan Hq values. For a fixed Hq, using Eq. (10) we 
identify intervals of /3 where cq is positive or negative. In 
these intervals, we find f3 values that satisfy Eq. (22) or 
Eq. (23) depending on the sign of cq. 


C. Interface model 


In Section fill B[ we have derived the JEM that treats 
the Kerr nonlinearity present in the core of the NSW 
in a simplified way [see Eq. @]- The used assumptions 
allowed us to obtain analytical formulas for the nonlinear 
dispersion relations and the field profiles of the nonlinear 
modes of the NSW. 

In this section, we will present the derivation of a 
model that is more numerical than the JEM but treats 
the Kerr-type nonlinearity in a more precise way. The 
permittivity of the nonlinear core in the frame of the in¬ 
terface model (IM) depends on all the components of the 
electric field and is described by 


£ 2 ( 2 :) = eiaix) + 02 [eI{x) + EI{x)] , (24) 


Moreover, there is no theoretical limitation of the values 
of the nonlinear permittivity change. 

In the IM, the solutions of Maxwell’s equations are 
sought numerically, as explained in the following. The 
field profiles inside the nonlinear core are found by nu¬ 
merical integration of Maxwell’s equations that couple 
the Ex and Ez field components. The novelty of our 
numerical method lays in the fact that, the parameter 
space where the solutions are being sought is reduced by 
a constraint that is expressed in an analytical form. Be¬ 
low the derivation of this constraint is presented and the 
numerical procedure of finding the nonlinear dispersion 
relations using the IM is shortly described. 


1. Analytical constraint 


The derivation of the IM starts from Maxwell’s equa¬ 
tions [Eqs. ([^]. In this approach the magnetic field is 
eliminated from these equations. The use of Eq. (2b) in 


Eqs. (2aI and (2cI gives [lUllTFlITfl 177] 




d{eEx) 

da: 


= PkoeEz. 


(25a) 

(25b) 


In the nonlinear core Eqs. 
formed to (for the detailec 


(|25a), (25b) and (24) are trans¬ 
derivation see Ref. [77]) 




-kl^{El + ElY + Co, 


(26) 


where the linear permittivity e/_ 2 , the nonlinear parame¬ 
ter 02 of the nonlinear core of the slot waveguide appear 
and Cq denotes the integration constant. 

In problems with a semi-infinite nonlinear medium 
the magnetic field Ely and its x-derivative vanish when 
± 00 . In such cases, these boundary conditions al¬ 


ee 

low us to set the integration constant Cq in Eq. (|26|) to 
zero [16] . Here we deal with a problem in which the non¬ 
linear medium is sandwiched between two linear (metal) 
layers, and therefore the nonlinear medium has a finite 
size. In this case, the integration constant can not be set 
automatically to zero. 


We compare the right-hand side of Eq. (261 with the 


square of the right-hand side of Eq. (25a I. The compar¬ 
ison gives 


{el/13^ - 262) El+ei,2 {eI + E^) + 
02/2 {El 


Et 


= Cq. 


(27) 


Equation (27) together with continuity conditions for the 


tangential components of the electromagnetic field will 
be helpful in finding the constraints reducing the param¬ 
eter space where the solutions of Maxwell’s equations are 
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sought, in order to find the dispersion curves for the NSW 
configuration. 

We use the continuity conditions for the tangential 
components of the electromagnetic field in order to re¬ 
late the values of the electric field components and 
Ez at the nonlinear interfaces to the values of the total 
electric field amplitude, defined as 

Ep EE ^El^ + El^, (28) 

where the additional subscript p S {0, d} denotes the x 
coordinate at which the quantity is calculated. The field 
distributions for the electric field components in the semi¬ 
infinite linear metal regions are found by solving linear 
wave equations for these components. The components 
of the electric field in the cladding metal layers are given 
by: 

1. In the left metal region (a; < 0 — layer 1 in Fig. [^: 


amplitude at this interface (as a parameter), the effective 
index f3 as an unknown and the opto-geometric material 
parameters which are known and fixed for a given NSW. 
Inserting Eqs. (31a) and (31b) into Eq. (27) taken at 
cc = O’*", we obtain the nonlinear dispersion relation at 
the left interface (x = 0 ): 


irr- 


2 e 


' 2,0 


(ei/3)^ 


(e2.o9i)^ + (ei/3)^ 


+ £ 1,2 + 


0^2 jz2 _ C'o 
2 ^0 " E 2 ■ 


(33) 


X = d , we obtain t 
interface {x = d): 


Inserting Eqs. (|32a|) and (|32b|) into Eq. taken at 
= d- 


re dispersion relation at the right 




0,0 


( 62 ,d®)^ + 


+ e /,2 + 


0.2 


2 " El 


(34) 


E,,i = (29a) 

Ez,i = (29b) 

2. In the right metal region {x > d — layer 3 in Fig.[^: 


Comparing the dispersion relations for the single in¬ 
terfaces given by Eqs. (33) and (34), we eliminate the 
integration constant Cq and obtain the final equation of 
the IM: 


Ez,o = 


-d) 

5 

-d) 


(30a) 

(30b) 


Only one exponential term is present in each of the ex¬ 
pressions so that the electric field decays exponentially 
for X -0 ±oo. 

Using the continuity conditions for the tangential field 
components {Ely and Ez) at both interfaces (cc = 0 and 
X = d), and Eqs. (2b) and (2c) in the linear layers, we 


can express the electric field components at the left in¬ 
terface {Ex,q, Ez,o) as a function of the total electric field 
amplitude at this interface Eq: 

Ex,o = {^iPEqI/[( e 2 ,oqil + (ei/3)^], (31a) 

EI,o = {£2 ,0^1 EqI/[{ e2,oqil + (ei/3)^], (31b) 

where e 2 ,o denotes the value of the nonlinear permit¬ 
tivity at the left interface of the core and is equal to 
e 2 U=o+ = e /,2 + Q: 2 ^'o- Similarly, we can express the elec¬ 
tric field components at the right interface {Ex,d, Ez,d) 
as a function of a total electric field amplitude at this 
interface E^: 

El, = {£sl3E,n/{e2,dq3f + {£313)% (32a) 

E^,d = {£2,dq3Ed)^/[{£2,dq3)^ + {£3(3)1, (32b) 


where e 2 ,d denotes the value of the nonlinear permittivity 
at the right interface of the core and is equal to £2 U=(/- = 
e/,2 


a 2 Ej. 


Equation (27) can now be rewritten on each interface in 


such a way that, it depends only on the total electric field 


(w) 


-2e 


2,0 


(=r)' 


- 2e 


2 ,d 


Eld? 

E2,ogiV+EidV 


+ e/,2 + %Eq 


E2,dq3? + E3d) 


2 + e /,2 + ?rE'( 


Ej 

El 


(35) 


Equation (35) represents the constraint that will be used 
in the numerical algorithm computing the dispersion 
curves presented in Section lTlI C 2[ in order to reduce the 
dimension of the parameter space where the solutions of 
Maxwell’s equations in NSW structures are sought. 


It is worth noticing that, in Eq. (35), neither the wave¬ 


length of light A nor the width of the waveguide core d 
appear. This means that the condition given by Eq. (35) 
is identical, regardless of the values of A and d. This 
condition depends only on the material parameters (ei, 
e/, 2 ) 02 , and £ 3 ) and the field intensities at both non¬ 
linear core interfaces {Eq and E,). The size of the core 
and the wavelength will appear in our next step — the 
numerical integration of Maxwell’s equations leading to 
the field profiles in the core of the waveguide and to the 
dispersion curves of the NSW. 

As stated before, Eq. (35) is not a dispersion relation 


for the slot waveguide modes but only a constraint that 
limits the parameter space where the solutions in the 
frame of the IM can be found. In order to obtain the dis¬ 
persion relations in the NSW, the field profiles in the core 
are found by numerical integration of Maxwell’s equa¬ 
tions. Maxwell’s equations are written as a set of cou¬ 
pled equations relating both electric field comp onent s E^ 
and Ez- These equations are derived from Eqs. (25a) and 
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(25b): 


iE, - 2a2E,El U - f) 

^ ea + 2a2E^ 


(36a) 

(36b) 


2. Numerical algorithm and nonlinear dispersion relations 

In this section, we present the description of the numer¬ 
ical algorithm that is used to find the dispersion relation 
for the modes of a NSW in the frame of the IM. This 
algorithm uses the shooting method | 1 S] to find solutions 
of Maxwell’s equations in the waveguide core and finally 
the nonlinear dispersion relation for the NSW. 

In a general case, the numerical procedure would be 
the following. First, for a given structure, we fix the 
parameters Eq, E^, and j3 and integrate Maxwell’s equa¬ 
tions [Eqs. (36)] in the waveguide core with the values 
Eq and j3 as initial parameters. The results of the in¬ 
tegration are the field profiles Ex{x) and Ez{x) inside 
the nonlinear core and, in particular, the computed total 
electric field amplitude at the interface x = d denoted 
by Next, we verify if the result of the numeri¬ 

cal integration fulfills the conditions resulting from the 
problem formulation: 


1. Is E 


(num) 


equal to the initially fixed value Ed^ 


2. Does the a;-derivative of the product of the per¬ 
mittivity and the transverse electric field e 2 Ex at 
the interface x = d~ have the correct sign? The 
condition for the correct sign reads 


E. 


x,d 


d{e2E^) 


dx 


— —koQsisBl, 


(37) 


and it is derived in Appendi x [b] and in Ref. [77]. 
The right-hand side of Eq. (37) is positive as ea 


of the metal is negative and all the other quanti¬ 
ties there are positive. Therefore, the sign of the 
derivative [d{e 2 Ex)/dx\\^^d- must be identical to 
the sign of E^^a- 

3. Do the components of electri c fiel d at x = d fulfill 
the conditions given by Eqs. (32a) and (32b)? 


Is these three conditions are fulfilled then the triplet {Eq, 
Ed, and /3) and the corresponding field profiles are ac¬ 
cepted as a genuine solution of our problem. 

In the general case described above, the parameter 
space where the solutions are sought is three-dimensional 
and it is spanned by Eq, Ed, and jd. However, we can 
separate the problem into two cases, where we will be 
able to simplify it and look for the solutions in only two- 
dimensional spaces. The two cases are : 

1. The case of symmetric {Ex,q = Ex,d) or antisym¬ 
metric solutions {Ex,q = —Ex,d) (for both sym¬ 
metric and antisymmetric solutions Eq = Ed) in 
a symmetric NSW (ei = £ 3 ). In this case, we look 
for the solutions of Maxwell’s equations in a two- 
dimensional space spanned by Eq and fd for each of 
the cases E„, q = EE^^d- For Eq = Ed in symmetric 
structures, Eq. ( |3^ represents an identity and it 
is satisfied for all values of (d. Therefore, Eq. (35) 


will not provide any help in further reducing the 
parameter space where the solutions are sought. 

2. The case of either the asymmetric solutions {Eq ^ 
Ed) in a symmetric NSW structure (ei = £ 3 ) or 
any solution in an asymmetric NSW (ci ^ £ 3 ). 
For these types of solutions, Eq. ( |^ is not an 
identity and results in a constraint on the three- 
dimensional space where the solutions are sought. 
Equation (35) is transformed to the form: 


P4;d +P 2/3 Tpo = 0, 


(38) 


where 


Pa — 2t2,d^\Ed{e2,o + ^ 1 ) ~ 2e2fie\El{e\ d + £ 3 ) + + ^i)i^2,d + £ 3)1 

_ 2 2 77’2 / 2 I 2 \ 2 2 77’2 / 2 i 2 \ i 

P2 — ^2,0^1^o{^2,d + ^ 3 ) ~ ^2,d^3^d\^2,0 + ^l) + 

2£2,o£2,d£i£3(£2.d£iF^o “ ^2,o^?,E‘ 1) - f[^\,Qti{e\^d + £3) + 4,d^3(£2.o ^ 

Po = ^2fi^2,d^l^3{f^sEd — SiEq + /), 

/ = £z, 2 (Ao + E^) + -e{Eq + E'^). 


£?)], 


(39a) 

(39b) 

(39c) 

(39d) 


This shows that Eq. (35) is satisfied only by a fi¬ 
nite set of (d values. Because we look for forward 
propagating nonlinear modes and material losses 


are neglected, the only physically meaningful solu¬ 
tions of Eq. (38) are the ones where fd is real and 
















(42) 


positive. Therefore, the two possible roots are: 


form: 


/ 3 ± = 


'-P2 ± - 4p4Po 

2P4 


(40) 


Only the real solutions among /3± are used further 
in the process of the resolution of the nonlinear 
problem. 


Using Eq. ( |35[ ), the three-dimensional space from 
the general case is now reduced to a two- 
dimensional space [one for each of the real effec¬ 


tive indices given by Eq. (40)] spanned by Eq and 
Ell- Therefore, instead of scanning a full three- 
dimensional space spanned by Eq and Ed, and /3 
we need to scan only a few two-dimensional spaces 
spanned by Eq and Ed corresponding to the phys¬ 
ically meaningful /3 values given by Eq. (401. In 


other words, for a pair {Eq, Ed) we just need to 
check if the held integration gives valid results (i.e., 
if conditions [T]|^ are fulhlled) for real values among 
/3±. If all the conditions are fulhlled then the triplet 
{Eq, Ed, and (3) and the corresponding held prohles 
are accepted as a genuine solution of our problem. 


IV. LIMITING CASES FOR TWO-LAYER 
STRUCTURES 

In this section we will derive the expressions for the 
dispersion relations in the limiting case of the single in¬ 
terface between a metal and a nonlinear dielectric. This 
will prove that our models reproduce already known re¬ 
sults for simpler structures. This limiting case disper¬ 
sion relations will also provide approximated analytical 
expressions for the propagation constant of highly asym¬ 
metric modes that resemble nonlinear plasmons on a sin¬ 
gle interface as shown in the following article |44j. 


A. Jacobi elliptic function based model 


In Sec. IIIB I we have stated that in case of a semi- 


inhnite nonlinear medium (a single interface between a 
metal and a nonlinear dielectric) the integration constant 
in Eq. Q should be set to zero as both the magnetic 
held Hy and its derivative tend to zero at inhnity. One 
way to hnd the dispersion relation for nonlinear waves 
propagating along a single interface is to use Eq. (ITOl. 


Setting Co = 0 in this equation we obtain an analytical 
formula for the effective indices 

(ei.2/ei)^g?-gi += 0. (41) 

Using the dehnitions of qk and 02 [see Eqs. § and 0 ] 
in Eq. (41), we hnd the analytical expression for the ef¬ 


fective index of a nonlinear wave at a single interface 
between a metal and a nonlinear dielectric in an explicit 


13 = 


\ 


eiei,2(e/,2 — ei) 

1 I 


^1,2 


2eoC€i^2 


Another way to hnd the dispersion relation for nonlinear 
waves propagating along a single interface is to set Cq at 
the later stage of the JEM derivation, namely in Eqs. (22) 
or (23) and the related parameters [Eqs. (15) or ( |17|), 
respectively]. In both cases this procedure yields [77] 
(using the expressions for the limiting values of Jacobi 
elliptic functions in the case of the parameter m = 1 
provided in Ref. |24j) 


tanh [fco 92 (d - a:o)] = . 

£392 


(43) 


Equation (43) describes the nonlinear dispersion relation 


for plasmon-solitons on a single interface (at x = d) be¬ 
tween a metal and a nonlinear dielectric. Equation (43) 


is equivalent to Eq. (42) from Ref. |T6|, which gives the 
dispersion relation for a single metal/nonlinear dielectric 
interface obtained using the held based model developed 
in Ref. jlhj . taking into account that: (i) in the frame 
of the JEM we used the assumption that £2 = e /,2 in the 
continuity conditions used to derive the nonlinear disper¬ 
sion relations (see Section III B 2 and (ii) in Ref. [TB] the 
interface is located at x = 0. 


Equations (42) and (43) give also approximated ex¬ 


pressions for the effective indices of highly asymmetric 
solutions in NSWs, as it will be proven by the numer¬ 
ical results presented in Ref. [44]. Highly asymmetric 
solutions are strongly localized on one of the interfaces 
and therefore the problem can be simplihed to a single¬ 
interface problem. A comparison of the approximated 


lutions of the JEM will 
Ref. mi. 


solution given by Eqs. (42) and (43) w ith the exact s o- 
3 e given in Sec. [section] [2] [ II] of 


It is worth noting that using Eq. (42) in the linear case 


( 2 ) 

{Hq —)■ 0 or 77,2 —>■ 0), we recover the dispersion relation 

for a linear surface plasmon propagating along a single 
interface {see Eq. (2.14) in Ref. |78] l 


P = \J eie/,2/(ei + e/,2)- 


B. Interface model 


(44) 


Equations (331 and (341 considered separately give the 


dispersion relation for a single interface between a metal 
and a nonlinear dielectric. In this case, the nonlinear 
medium is semi-infinite which means that Cq must be 


set to zero (see discussion in Sec. Ill C 1). Setting Cq = 0 
in Eq. (331 yields 






(e2,o<?i)^ + (ei/3)2 


+ e /,2 + ^E^ = 0 . 


(45) 
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This equation can be solved analytically for /?. The solu¬ 
tion depends on the parameters of the structure (ei, ei^ 2 , 
a 2 , £ 3 ) and the electric field amplitude at the interface 
Eq and is given by 


/? = 


<=1<^2,o(^1,2 - Cl + 

(^ 2,0 + “ ‘^^ 1 ^ 2,0 


(46) 
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Only the positive root is considered because we are in¬ 
terested here in forward propagating waves only. Equa¬ 
tion (46) can be compared to Eq. (11) in Ref. [74] and 
to Eq^l4) in Ref. [10] derived for the case of a single 
metal/nonlinear dielectric interface. 

The solution for a single-interface problem provides a 
good approximation (in terms of effective index /3 and 
the field profiles) for first-order highly asymmetric modes 
in the NSWs whose field profiles are mostly localized 
on one interface of the core (see also the discussion in 
Sec. IV AI. These solutions are invariant with respect to 


the waveguide width as they interact strongly only with 
one of the core interfaces. More comments and illustra¬ 


tions of this property will be presented in Section [sub- 
2 2 lIB]in Ref. [44], where we discuss the results 


section 


for the symmetric NSWs. In the linear limit a 2 E^ 
Eq. (]4^ transforms to Eq. (44), as expected. 


Appendix A: JEM 

In this appendix we present the derivation of the so¬ 
lution of Eq. ([ 8 ) that leads to the expressions for the 
field profiles in the nonlinear core in the frame of the Ja¬ 
cobi elliptic function based model (see Sec. jlllBl] ). The 
derivation is presented in separated cases depending on 
the signs of the integration constant cq and of the pa¬ 
rameter. This separation allows us to work with positive 
quantities only in each of the cases, which facilitates the 
choice of ambiguous signs. 

In this appendix we first divide the problem according 
to the sign on For each of the possible signs we di¬ 
vide the problem according to the cq sign. This choice 
is dictated by the fact that in the case of negative we 
have to consider only the case of positive Cq as it can be 
seen from Eq. (10). 


V. CONCLUSIONS 


1. The case gf > 0 


a. The subcase co > 0 


We have presented two complementary models based 
on Maxwell’s equations to study the properties of the sta¬ 
tionary TM solutions in planar nonlinear structures made 
of a focusing Kerr nonlinear dielectric core surrounded 
by semi-infinite metal regions. The first model uses a 
simplified treatment of nonlinearity that takes into ac¬ 
count only the transverse component of the electric field 
in the nonlinear term and assumes that the nonlinear 
refractive index change is small compared to the linear 
refractive index. It provides an analytical description of 
both the field profiles in the whole waveguide using Ja¬ 
cobi elliptic special functions and the nonlinear disper¬ 
sion relations. These features allow one to study rapidly 
and accurately the properties of the NSW structures as a 
function of their opto-geometric parameters. The second 
model takes into account the full dependency of the Kerr 
nonlinear term on all electric field components and no 
assumption is required on the amplitude of the nonlin¬ 
ear term. It allows us both to determine the domain of 
validity of the first model and to investigate with accu¬ 
racy the effects of high nonlinearities. These two models 
prove their usefulness in the next article (Ref. [U]) where 
we provide a complete and accurate study of nonlinear 
slot waveguides. It is worth mentionning that these mod¬ 
els can be extented to more complicated geometries and 
refractive index forms. 


At first, we consider the case where > 0 and Cg > 0 
and find the solutions of the nonlinear wave equation 
[Eq. ([^] for this case. Equation (|^ was transformed 
into Eq. (12) and its left-hand side can be expressed in the 


(Al) 


form of the integrand of an elliptic integral (see Ref. 

dH2 _ + , /l H.. 

where the parameters 7 and S were introduced. In order 
to relate the newly introduced parameters with parame¬ 
ters A, Q, and cg, we compare the expressions unde r th e 
square-root on the left-hand sides of Eqs. (12) and ( [Al] ). 
This comparison results in 

72 = (^A2Q2+4Aco - AQ)/2, (A2a) 

= (^A2Q2+4Aco + AQ)/2, (A2b) 


where the choice of the sign in front of the AQ term is 
dictated by our assumption that the magnetic field Tfy 
is real [see the assumption made in Eq. ([^ and Ref. [77] 

for more details]. _ 

Integrating Eq. (]A1[), we obtain 


Jh2{0) 1/(72 -ki7|)((52 - Hi) 



(A3) 
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The integral on the left-hand side of Eq. (A3) can be 
separated into two integrals: 


Using this fact, Eq. (A31 yields 


rH2{x) 

'^ 2 ( 0 ) 


lH2iO) JH2 (x) 


(A4) 


AH 2 




/ff2(o) - Hi) Jh2{x) \/(7^ - Hi) 



(A5) 


r 


Multiplying Eq. (^A5| by -|- <5^ and using formula 
17.4.52 from Ref. |24) . gives 


cn 


-1 


H2{0) 


— cn 


-1 


H2ix) 


m 


= ± 


'jx, (A 6 ) 


where cn ^ (m|to) is the inverse of the Jacobi elliptic func¬ 
tion cn(M|TO) and 


m = 

s = 7 ^ -I- S^. 


(A7a) 

(A7b) 


Jacobi elliptic functions are defined with the argument u 
and the parameter m [24) . Reorganizing the terms and 
appl ying the Jacobi elliptic function cn to both sides of 
Eq. (A 6 ), results in the expression for the magnetic field 


in the core of a NSW in the case where both cq and ql 
are positive: 

H 2 {x) = Jcn|=F \Js/A [x — Xf:))\m^ , (AS) 

xq = —sj A/scn“^ [iJ 2 ( 0 )/J|m] . (A9) 

The uncertainty about the sign in front of the square- 
root in the argument of the cn function in Eq. (AS) can be 


resolved by the analysis of Eq. (A3) using the properties 


of the Jacobi elliptic function cn and the field continuity 
conditions at the interface x = 0. Without loss of gener¬ 
ality of the obtained results, we assume that in iJo > 0 
and d > 0. In this case, the proper choice of the sign in 
front of the square-root in Eq. (AS) is the bottom plus 
sign (see Ref. m for more details). Finally, the field 
profile in the nonlinear core for the subcase La is given 
by 

H2{x) = Jcn I a/s/A {x — xo) |to| . (AlO) 


b. The subcase co < 0 

Here we consider the case where <72 > 0 and Cq < 0. In 
order to work with positive quantities only, we substitute 
Co by — |co| in Eq. (12) and obtain 


diJa 


\/—A|co| -l- AQHl — Hi 



The left-hand side of Eq. (All I is then expressed in the 


form of the integrand of an elliptic integral (see Ref. [24] 1: 


AH 2 

v/(7^ - Hl){Hl - J2) V A 


(A12) 


Comparing the expressions under the square-root on 
the left-hand sides of Eqs. (All) and (A12) allows us to 


find the relations between the new parameters 7 and 5 
and the parameters A, Q, and Cg: 


7 ' = (AQ -f ,jA^Q^-AA\cc>\)/2, 
= [AQ - VA2Q2_4yi|co|)/2. 


(Al3a) 

(A13b) 


The choice of the signs in front of the square-roots is 
Eqs. (A13| is arbitrary due to the symmetric role of 7 


and 5 in Eq. (A12 1 and due to the fact that = AQ — ^'^ 


The expression that appears under the square-root in 
Eqs. (A13) is a difference of two positive quantities. In 


order for 7 and J to be real (which ensures that the mag¬ 
netic field H 2 is real), the quantity under the square-root 
must be positive or equal to zero. Writing expression un¬ 
der the square-root explicitly using the definitions of A, 
Q [Eqs. ([1§], and Cg [Eq. ( |l0| )] we obtain 


A2g2 4^cg = 

(<72 ~ ®2JJg) + 2a2iLg (e2|a;=g+/ei) ql, (A14) 

which is greater or equal to zero because both terms in 
the sum are greater or equal to zero. This proves that 7 
and 6 are real quantities. 

Having found the expressions for 7 and 6, we proceed 
with the derivation process. Integrating Eq. (A12) gives 


_dH2__ [T ^ 

Ih2{o) ^( 7 ^ - Hl){Hl - (52) V A 


(AI5) 


This time, the integral on the left-hand side of Eq. (A15) 
is separated in the following way: 


rH2{x) 


’H2(0) 


lH2iO) JH2 {x 


(All) Inserting Eq. (A16) into Eq. (A15), we obtain 


(A16) 
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dm 


dm 


Ih2{0) V(7^ - - S^) Jh2(x) a/(7^ - - 6^) 



(A17) 


Multiplying Eq. ( |A17 ) by 7 and using formula 17.4.44 
from Ref. [21] yields 


dn 


-1 


m{o) 


7 


m 


— dn 


-1 


m{x) 


7 


m 



(A18) 

where dn ^ is the inverse of the Jacobi elliptic function 
dn and 


= (7^-^W- 


(A19) 


Reorganizing the terms and applying the Jacobi elliptic 
function dn to both sides of Eq. (AI 8 ) results in the ex¬ 


pression for the magnetic field in the core of a NSW for 
positive and negative cq: 

772 ( 2 ;) = 7 dn I (a; - Xq) Im| , (A20) 

xo = -\J dn“^ [i72(0)/7|m] . (A21) 

where the sign in front of the square root is chosen in a 
similar manner as in the case Cq > 0 . 

2. The case gl < 0 

H ere we consider the case of < 0. As stated in 
Sec. IIIB below Eq. (10), in this case, the integration 


constant cq takes only positive values. In the following, 
we find solutions of the nonlinear wave equation [Eq. ([^] 
for this case. 

In or der to work only with positive quantities in 
Eq. (Ill, in the case of negative q^, we will substitute 
by its negative absolute value — l^ll- This substitu¬ 
tion transforms Eq. (11) into 


dm 


y^co - fc2|g||i7| - k§fHi 


= ± dec. 


(A22) 


We redefine Q given by Eq. (13a) to be positive. In the 
case II Q is defined by 


Q = ko\ql\ 


(A23) 


and A is still defined by Eq. (13b). Using this definition, 


Eq. (A22| can be written in the form: 


= ± dx. 

v/Aco - AQHl - Hi V A 


(A24) 


We rewrite Eq. (|A24|) in the form of the elliptic integral: 

(A25) 


dm /I 

V(72+i7|)(<52-Jj|) V A 


Comparing the expressions under the square-root on left- 
hand sides of Eqs. ( |A24 1 and (A251 allows us to find the 
relations between the new parameters 7 and 6 and the 
parameters A, Q, and cq: 

7 ^ = (v/A2Q2+4^co + AQ)/2, (A26a) 

= (y/A2Q2+4Aco - AQ)/2. (A26b) 

The following derivation is exactly the same as in the 
case La but we have to keep in mind that the definitions 
of 7 and S are reversed [compare Eqs. (A2) and (A26l] 


and that Q is defined now by Eq. (|A23|). Equation (|A25|) 
is integrated and yields a formu 
Eq. ([Al: 


a that is identical to 


rH2(x) 


dm 


lH2iO) y /(72 + 7L2)(^2_^2) 



(A27) 


By an analogy to the case presented in Sec. [A 1 a[ the 
final expression for the field profile in the nonlinear core 
for < 0 has the form: 

m{x) = (5 cn I -\/s/A (x — xo)|m.|, (A28) 

Xq = —\/A/scn“^ \^m{d)/5\m\ . (A29) 


3. 


Summary and unification of the expressions for 
field profiles in the nonlinear core 


In Sec. [BIB Ij we have split the problem into two cases 
depending on the sign of the integration constant cq in¬ 
stead of three cases presented in this appendix. Below we 
show how to obtain the splitting presented in Sec. BIB Ij 
from the results derived in Appendices [A Ij an d |A^ 

We notice that the case presented in Secs. [A 1 aj and 
A 2 can be merged into one case. The sole difference 


between the expressions obtained in the case [A 1 aj and 
A 2 is in the formulas for the parameters 7 and 5 [com¬ 
pare Eqs.(A2) and Eqs.(A26)]. The sign in front of Q 
is reversed. This difference is caused by our choice of 
the definition of the Q parameter [compare Eqs (13a) 
and (A23)]. This choice allowed us to work with positive 
quantities only which simplified the procedure of deter¬ 
mination of the ambiguous signs. During the derivation 
of the nonlinear dispersion relations we will not encounter 
such problems any more and we can work with Q that 
is either positive or negative. Therefore, without loss of 


generality we unify the cases presented in Secs. Ala and 
|A 2| to a common case where cq > 0. Now the distinc¬ 
tion between the cases will be based on the sign of the 
integration constant cq. 
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In the derivation of the dispersion relations presented 


in Sec. IIIB2 we have considered two cases: I (cg > 0) 


for which we will use Eqs. (13), (|A2 


(AlO) [presented as Eqs. (14) and (|15) in Sec. 


used. 


1 , ( 
p5) 

A7 

), (AC 

), and 

in Sec. I 

IlM]; 

(AH 

5), (A19)- 

-(A2I) 

Sec. 

III B I will be 


In the nonlinear medium (layer 2) Eq. (B2) is written in 
the form 


A, = 


1 

kol3e2 dx 


(B4) 


Inserting Eqs. (B3) and (B4) into Eq. (Bl) yields 


Appendix B: Derivation for the analytical condition 
in the numerical algorithm 

In this appendix we present the derivation of the con¬ 
dition in Sec. El C 2| used in the numerical algorithm 
of the IM. We start the derivation from the continuity 
condition for the ^-component of the electric held at the 
interface x = d: 


1 d{e2E^) 

e2,d dx 


dE^ 


dx 


X — d+ 


(B5) 


Next we use the continuity condition for E[y at the inter¬ 
face X = d 


Hy^2\x = d~ Ely^'i\x = d+^ 


(B6) 


which with the help of Eqs. (2b) and (30a) becomes 


E 


z,2\x—^ 


=d- — E. 


'z,Z\x—d'^ • 


(Bl) 


i2,dEx,d = 


(B7) 


The longitudinal held component E^ can be expressed 
with the help of Eq. (25b) as 


A, = 


1 d{eEx) 
kof]e dx 


(B2) 


where we consider isotropic medium for which = Cz = 
e. In the case of the uniform linear medium (as it is the 
case for layer 3), Eq. (B2) simplihes to 


Ez = 


1 dEx 
ko/3 dx 


(B3) 


Equation (B5) is now multiplied by Eq. (B7) to give 
d{e2Ex) „ dE, 


E. 


x,d 


dx 


— ^sBx 


C=d~ 


dx 


(B8) 


c—d+ 


Computing the derivative on the right-hand side of 
Eq. (B8) using Eq. (30a) gives [see Eq. (37)] 


Ex 


d{e2Ex) 


dx 


— —koqsesBx- 


(B9) 
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